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I.  STABILITY  AND  CONTRA CTTVITY  OF  MULTISTEP  AND  ONE-LEG  METHODS 
A.  Background 

With  a  linear  multistep  formula, 

y-0  7-0 


normalized  by 


k 


7-0 


(2) 


we  can  associate  two  integration  methods  for  systems  of  ordinary  differential  equa¬ 
tions 


i  -  /«,*) 


(3) 


namely  the  linear  multistep  (MS)  method 


k  k 

2  aj**+j  ~h2  0/('n+7’W  ■  0 

7—0  7—0 


(4) 


and  its  one-leg  (OL)  "twin" 


2  VW  -  W  2  */W  2  -  °- 

7—0  7—0  7—0 


(5) 


The  formula  (1)  is  said  to  be  stable  at  q  ■  h\  if  all  solutions  {*„}  of  the  difference 
equation,  generated  by  applying  (1)  to  the  test  equation 


i  Ax,  A  complex,  constant. 


(6) 


(for  which  the  MS-  and  OL-methods  are  identical)  remain  bounded  for  a  given  A>0  as 


-  1  - 


This  is  the  case  iff  the  characteristic  polynomial 

X«)  -  p(f)-*«(f)  (7) 

satisfies  the  "root  condition"  (i.e.,  ail  roots  f,  of  x(f),  <  ■»  1  satisfy  |  f,  |  $  1  and 
If, |  »  l^f,  is  simple).  Here 

k  k 

p(f )  -  s «/.  «(f )  -  2  0/-  (») 

y-0  y-0 

The  stability  region  5  is  the  set  of  all  q  at  which  the  formula  is  stable.  The  formula  is 
said  to  be  A-,  A0-,  and  A^-stable,  if  0_:  «  {q  |  Re  q<0},  R_:  ■  (qlq<0),  or 
lq  -  <*>}  are  contained  in  S,  respectively. 

The4ormula  (1)  is  said  to  be  contractive  [8]  at  q  with  respect  to  a  given  norm 
|*|  if  |X„_,  |  >  ||.Y/I|  for  all  n>k,  where  XH:  -  (xn.k+\yXn_k+2,..  •.*„)•  The 
contractivity  region  K  (with  respect  to  |  •  J  )  is  the  set  of  all  q  at  which  the  formula  is 
contractive.  The  formula  is  said  to  be  A-,  AQ-,  or  A^-contractive,  if  0_,  R_,  or 
{q  m  oo}  »re  contained  in  K,  respectively.  Clearly,  contractivity  at  q,  A-contractivity, 
etc.,  imply  the  corresponding  stability  properties,  since  KiS.  Contractivity  is  a  useful 
tool  for  obtaining  stability  results  for  nonlinear  problems  or  equations  with  variable 
coefficients,  and  for  variable  integration  steps,  as  well  as  for  identifying  novel,  very 
robust  integration  formulas.  A  survey  of  contractivity  results  is  given  in  [17]. 

B.  Contractivity  Results 

1.  A0-  and  A(a)-contractivity  results. 

The  formula  (1)  is  A^-contractive  with  respect  to  the  f^-norm  iff 

*-l 

y:  -  0k-2  I  P.  I  >0.  Subject  to  the  normalization  (2),  one  has  y$l.  We  proved  [8] 

y-0 

that  for  any  p  there  exist  formulas  of  order  p  which  are  A(a)-contractive  with  respect 
to  the  max  norm  for  any  a  smaller  than,  but  arbitrarily  close  to  */2  and  whose 
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measure  of  A^-contractivity,  y,  is  arbitrarily  close  to  unity.  (Note  that  since  A- 
contractivity  implies  A-stability,  this  result  for  pfi  is  as  strong  as  it  can  be).  In 
general,  however,  the  length  k  of  these  methods  grows  approximately  like  3p. 

We  derived  [14]  the  two-parameter  class  of  all  two-step  (k  ■  2)  second-order 
(pm  2)  formulas  which  are  A0-contractive  w.r.  to  the  maximum  norm  1*1,  for 
arbitrary  step  ratios  r  «  rH  m  hH/hn_ This  generalized  the  corresponding  results  for 
uniform  steps  given  in  [8].  For  the  test  problem  x  m  \(t)x  with  arbitrary  A(/)<0,  the 
OL-implementation  of  any  variable-step  A0-contractive  formula  provides  an  A0-stable 
method  for  arbitrary  step  sequences  [/»„}. 

2.  A-contractivity  for  two-step  methods 

We  proved  [8]  that  a  MS-formula  which  is  A0-contractive  with  respect  to  the 
max  norm  is  A-contractive  if  and  only  if  it  satisfies  the  constraint 

*2  K«J  +  /*>)/(«*  +  pfa)]l/2Sh  0<n<  +  »■  (9) 

i- o 

In  general,  (9)  is  not  easy  to  analyze  but  for  the  2-parameter  family  of  all  MS- 
formulas  with  p>2,  defined  by 


1(— 1  +  h,),  a,—  — hj,  «2 ^(1  +  hj). 


(10) 


+  fix 


1(1 -h0),  02-l(l  +  bo  +  h,), 


the  condition  (9)  is  satisfied  iff  [8] 

h0  -  1  -  b\ ,  Osh, Si. 


(ID 
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As  stated  in  [8],  the  special  formula 


•  0 


(12) 


minimizes  a  bound  for  the  global  truncation  error  derived  in  [9],  over  all  A-contractive 
p  m  k  m  2-formulas.  (See  subsection  I.  D.l  hereafter.) 

We  generalized  the  one  parameter  class  of  all  p  «  k  ■»  2-formulas  which  are 
A-contractive  w.r.  to  1*1.  to  arbitrary  step  ratios  r  »  rn  —  hH/hn_x  [14].  It  is 
defined  by 


— — (1  -  v),  0O  «-^-  +  ^-  — 


1  +r  2(1  +r)  2(1  +/•) 


■u-U-K l-v)],  +  W 


(13) 


x2,n 


1--~7(1~V). 


1  +  2r~  1  v  _  r  t.2 

l+r  2(l+r)  ~  2(l+r)  * 


where  the  formula  is  to  be  read  as  X  a, 
ter  v  satisfies 


\^0 ajJixn -2+j~h„^Pj'„X„-2+j 


0,  and  the  parame 


^—Sv^l.  (14) 

We  proved  [14]  that  any  formula  with  p  m  k  «  2  is  A-contractive  w.r.  to 
I  •  |  m  iff  it  is  G-stable,  i.e.,  A-contractive  w.r.  to  some  G-norm  (a  norm  associated 
with  a  positive  definite  quadratic  vector  form).  For  any  given  A-contractive 
p  m  k  m  2  formula  with  r  ■«  1  (uniform  steps),  specified  by  a  parameter  value  v  -  v,, 
Osvj^l,  A-contractive  extensions  to  non-uniform  steps  (r#l)  can  be  defined  in  such 
a  way  as  to  keep  the  G-norm  fixed  w.r.  to  n  [14).  If  this  is  done  then,  for  an  arbitrary 
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step  sequence  {*„},  G-stability  is  assured  for  any  dissipative  (monotone  negative) 
nonlinear  system. 

We  proved  [13,14]  that,  except  for  borderline  cases,  the  A-stable  and  A- 
contractive  p  m  k  —  2-formulas  can  be  identified  by  a  local  analysis  near  q  m  0  and 
for  small  non-uniformities  of  the  grid.  More  specifically,  for  uniform  and  nonuniform 
grids,  stability  in  any  arbitrarily  small  left  half-neighborhood  of  q  —  0  is,  roughly 
speaking  equivalent  to  A-stability  (for  the  latter  in  the  formal  algebraic  sense  that 
p(f)  -  qo(S)  satisfy  the  "root  condition"  (r.c.)  for  all  q.  Re  q<0).  Furthermore, 
requiring  that  the  r.c.  remain  satisfied  for  r  —  r„  —  1  +  *  under  arbitrary  perturbations 
around  oO  singles  out  all  A-contractive  methods  with  uniform  steps. 

3.  A-contractivity  for  formulas  of  arbitrary  length 

In  [16]  we  generalized  the  derivation  of  the  variable  step  A-contractive 
second-order-formulas  mentioned  above  from  km  2  to  arbitrary  length.  First  of  all, 
the  formula  (1)  is  A0-contractive  iff  [8] 

-*j>0,jmO,...,k-  1, 

(15) 

0*-2  10,1*0, 

J- 0 


and  at>0  (the  latter  follows  from  the  first  of  relations  (13)  by  consistency  and  the 
assumption  that  «k#0).  Subject  to  (15)  and  for  a^ytO,  the  formula  is  A-contractive  iff 


,-0 


where  q  -  iy,  y  is  real  and  ij  —  y1.  By  consistency,  F(n)  ■  —  i\F  +  0(n2)  for  n-*0, 

2 


and 


F-  -  £  (pj/*j)ZO 

7-0 


07) 


Is  thus  necessary  for  (16)  to  hold.  We  proved  a)  that  the  only  p  <-  2-formulas 
satisfying  (17)  are  defined  by 


0;  “  “  o*1 . fc 


(18) 


where  0,  »  0yfl:  «  and  A2:  -  ~-2$l__jay,  and  b)  that  (17)  is  also 

sufficient  for  (16)  to  hold  and  thus  (15)  and  (17)  together  define  the  set  of  all 
A-contractive  formulas.  Consistency  requires  that 


k—2 


-1-2  («*_,- !)(-«,). 

7-0 

*fc-l  *  “  ^1~  2  **-/“*>) 


(19) 


and  thus  the  A-contractive  formulas  are  in  one-to-one  correspondence  with  the  points 
of  a  simplex,  in  the  space  spanned  by  the  ( k  —  1)  free  parameters  a0>...,ak_1,  whose 

vertices  are  the  origin  and  f-ay  m  1  /9k  ■,  ai  -  0,  i+j;  j  -  0 . k  -  2).  These 

vertices  represent  respectively,  as  degenerate  y'-step  formulas,  the  Trapezoidal  Rules 
with  step  lengths  t„  —  tn_j,  j  —  0 .....Hr. 

To  illustrate  the  significance  of  the  above  theorem  we  remark  that  stability  can 
be  guaranteed  for  arbitrary  step  sequences  if  at  each  integration  step  any  one  of  the 
A-contractive  OL-methods  is  applied  to  x  ■  X(x)x  with  any  A(f),  Re  X(f)£0.  This 
behavior  is  in  marked  contrast  to  that  of  the  popular  two-step  second-order  backward 
differentiation  formula  which  can  be  destabilized  by  certain  combinations  of  variable 
steps  and  coefficients,  as  exemplified  hereafter: 

-6  - 


For  exponentially  varying  grids,  i.e.  n  m  and  for 

\(r)  m  rA,X,[/i  +  (r  —  l)f]  *,  the  formula  and  the  difference  equation  generated  by 
applying  a  one-leg  method  to  i  «  X(r)x  have  constant  coefficients  and  stability  can  be 
analyzed  algebraically  [14].  For  this  problem  with  increasing  steps  (r>l),  the 
backward  differentiation  formula  (BDF)  with  p  -  k  m  2  loses  A-stability  [12  —  14]. 
For  r<l  the  BDF  is  A-stable  but  when  applied  to  purely  oscillatory  problems 
(Re  X|  ■  0)  it  is  too  strongly  damped.  The  A-contractive  p  -  k  ■»  2  one-leg  me¬ 
thods,  on  the  other  hand,  are  A-stable  for  this  problem  for  any  r  and  accurately  reflect 
the  constant  amplitude  of  purely  oscillatory  solutions  [12,13]. 

As  a  generalization  of  the  corresponding  result  for  k  <■  2  mentioned  above  we 
proved  in  [16],  for  arbitrary  k,  that  the  A-contractive  formulas  for  uniform  steps  can 
also  be  derived  without  invoking  the  contractivity  concept,  but  instead  by  analyzing  a 
necessary  algebraic  condition  for  A-stability  for  mildly  non-uniform  grids.  The 
condition  in  question  is 


lf,(*)l<l  (20) 

for  q  »  iy,  y  real,  |_y|  sufficiently  small,  where  tv(q)  denotes  the  principle  root  of  the 
characteristic  polynomial  2(o j-qPj)tJ  (the  one  defined  by  f,(0)  «  I).  The  coeffi¬ 
cients  0y  are  associated  with  a  mildly  variable  grid,  i.e.,  one  for  which 
Oj  m  j  +  tj,  |  tj  |  <<l.  By  a  perturbation  analysis  with  respect  to  the  tj  one  can  show 
that  the  only  p  *  2-formulas  for  constant  steps  (tj  -  0)  which  can  smoothly  be 
extended  to  arbitrary  tj,  ( tj  |  small,  in  such  a  way  that  (20)  remains  satisfied  to 
first-order  in  the  e;  are  the  A-contractive  formulas. 

C.  Nonlinear  Stability  Results  Based  on  Contractivity 

In  [9]  we  gave  a  boundedness  result  for  solutions,  generated  by  an  A(0)- 
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contractive  formula,  of  a  system  of  the  form 


x  -  A(fpc)x  +  d(tpc). 


(21) 


where  A  m  (ay(/,x))  is  a  real  matrix  with  au< 0  and  a  certain  amount  of  diagonal 

dominance  and  d(t*)  is  bounded.  Subsequently,  in  (13],  we  gave  a  sharper  result  and 

a  more  transparent  proof.  Specifically,  we  showed  that  if  we  apply  to  (21)  a  one-leg 

method  which  is  both  contractive  at  q  -  0  and  (strictly)  contractive  at  q  »■  »  with 

*-i 


respect  to  the  /^-norm  (i.e.  if  oy< 0,  j  ■»  0 -  \,ak>0  and  y:  —  /)*— 2  \  Pj  I  >0, 

respectively)  then  we  have  input-output  stability  in  the  f^-norm  provided  B  —  (btJ)  is 

k 

a  M-matrix;  here  -  -  a,,y  and  bn  *  -(2  I  p.  |  sup  |  a.. (at, ox)  | .  The  two  results 

1  j- 0  1  trX  1 

mentioned  in  this  paragraph  are  the  only  nonlinear  stability  results  we  are  aware  of 
which  are  based  on  contractivity  analysis  and  which  do  not  assume  monotonicity  of 
the  nonlinear  system  but  rather  a  type  of  diagonal  dominance  of  A.  Note  that  this 
result  is  of  a  novel  type  in  the  sense  that  it  is  valid  for  all  h>0  without  requiring 
A-stability.  Earlier  nonlinear  results  based  on  stability  and  relevant  for  stiff  problems 
either  required  A-stability  (and  thus  apply  only  to  methods  with  p<2)  or  they  were 
valid  only  for  sufficiently  large  h. 


D.  Special  Methods  and  Aspects  of  Implementation 

1.  Special  contractive  methods. 

A  special  A0-contractive  formula  with  p  •  k  -  2  was  selected  in  [9]  which 
minimizes  a  bounded  of  the  global  error  produced  by  applying  any  A(0)-contractive 
OL-method  to  the  model  problem  x  -  \(t)<-a< 0,  />0.  The  formula  in 

question. 


xii+2-x«+l 


•  o. 


(22) 
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is  of  Adams  type.  We  proved  [9]  that  the  one-leg  implementation  of  this  method 
remains  A0-contractive  (and  thus  A0-stable)  for  any  variable  step  sequence  whatsoever 
as  well  as  for  r-dependent  X.  For  r  ■  rH  -  hn/hn_i,  its  coefficients  are  defined  by 

®0,n*0’  al.n*  ~  *»  a2,nml’  ^O.n  “  r/<2  +  2r>*  01,„*0*  and  h.n  *  (2  +  O/2  +  2 r). 

By  minimizing  the  same  objective  function  over  the  class  of  all  A-contractive 
p  m  k  m  2-formulas  with  uniform  steps  we  found  [8,14]  the  special  A-contractive 
formula  given  by  (12)  which  is  associated  with  v  ■  v,  ■  -j.  This  formula  was 
generalized  [14]  by  extending  Us  minimality  property  to  the  variable  step  case.  For 
arbitrary  r  it  i»  defined  by  v  «  2r/(2r  +1)  and  has  coefficients 

r2  R  r(l+r) 

»  Pq  ■“ 


(l+r)(l+2r) 


(l+2r) 


H-r 

l+2r 


Pi- *0 


(23) 


\+2r+2r 
(1  +r)(l  +2r) ' 


l+2r+2r* 
(1+2  r)2 


2.  Formula  selection  strategy  for  global  error  control. 

Given  any  sequence  of  steps  {An}  whose  ratios  rn  •«  hjhn_x  are  bounded 
away  from  zero  and  infinity,  and  any  X(r)  satisfying  Re  \(t)<  -  a,o>0,  we  derived 
[13]  a  simple  strategy  (local  in  t)  for  selecting  a  sequence  of  A-contractive  one-leg 
methods  with  p  ■  k  m  2  in  such  a  way  as  to  minimize  an  upper  bound  of  the  global 
truncation  error  produced  by  integrating  x  —  \(t)x. 
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3. 


Multistep  and  one-leg  methods  for  mixed  differential  algebraic  systems 
Systems  of  the  form 


FUJcj tjO-O 
GUjc  j)  m  0 


(24) 


are  often  encountered  in  applications.  Here  F  and  x  are  vectors  of  the  same  dimen¬ 
sion  and  so  are  G  and  y.  In  [10]  we  proposed  and  analyzed  one  MS-implementation 
and  two  OL-implementations  of  the  formula  (1).  The  particular  methods  of  this  type 
associated  with  the  variable  step  version  of  the  A(a)-contractive  Adams  formula 
defined  by  (22)  was  shown  to  be  second-order  accurate,  both  in  the  "differentiated 
variables"  x  and  in  the  "state  variables"  y.  All  three  implementations  were  tested 
numerically  on  a  nonlinear  diode  circuit  problem.  The  numerical  results  were  in  very 
good  agreement  with  the  theoretical  findings. 
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n.  NONLINEAR  STABILITY  OF  HIGH  CRDER  METHODS  FOR  STIFF 

EQUATIONS 

A.  Background 

High-order  linear  multistep  methods  (LMM)  have  computational  advantages  in 
the  numerical  integration  of  stiff  problems  whose  solutions  are  sufficiently  smooth. 
One  usually  assumes  in  such  situations  that  the  spectrum  of  the  Jacobian  of  the 
nonlinearity  stays  within  a  sectorial  set  and,  correspondingly,  one  uses  high-order 
A(a)-stable  schemes.  The  techniques  used  in  studying  the  nonlinear  stability  of 
low-order  methods,  e.g..  A-  and  G-contractivity  are  not  applicable  to  the  high-order 
ones  and  a  corresponding  nonlinear  stability  theory  has  been  lacking.  We  have 
introduced  [IS]  a  new  technique  to  investigate  the  behavior  of  the  numerical  errors  in 
high-order  LMM  when  applied  to  parabolic-like  nonlinear  problems.  We  consider  the 
global  error  equation 


ye„  +  h[f(xn  +  <?„)-/(*„)]  -  hqn,  (25) 

where  y  is  the  /,-sequence  defined  by  po~l({)  m  X 1  xn  is  the  approximate 

j 

solution,  qn  is  the  local  error,  en  the  global  error,  h  the  uniform  time  step,  /(•)  the 
nonlinearity  in  the  stiff  differential  system  and  •  denotes  convolution.  Instead  of 
scalar  multiplication  of  (25)  with  e„,  we  multiply  by  M*%>  where  p  is  an  appropriately 
chosen  sequence  called  the  "multiplier"  and  is  defined  hereafter.  This  procedure 
yields  two  "energy-like"  error  terms,  namely  a  quadratic  term  describing  the  interac¬ 
tion  between  the  multiplier  and  the  method,  and  a  more  general  nonlinear  term 
combining  the  problem  with  the  multiplier.  To  control  the  global  error,  both  terms 
have  to  be  positive  and  our  stability  theory  finds  conditions  for  such  positivity.  The 
results  naturally  subdivide  to  cover  three  areas:  the  relationship  between  the  multistep 
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method  and  the  multiplier,  the  interaction  of  the  nonlinearity  with  the  multiplier,  and 


I* 


Since  it  is  easy  to  see  that  (27)  implies  A(a)-stability,  Theorem  1  may  be 
interpreted  as  saying  that  the  linear  stability  of  an  A(a)-stable  method  can  be  "seen" 
through  the  multiplier. 

2)  A  quantitative  "uncertainty  principle"  was  derived  showing  the  incompati¬ 
bility  of  extreme  stability  and  accuracy  of  a  method.  The  error  in  a  method  may  be 
measured  by  the  size  of  the  Fourier  transform  of  the  associated  Peano  kernel  K,  and 
its  stability  by  y~l  where  (1,— ij,  0, ...)  is  the  special  Popov  multiplier  associated  with 
the  method.  In  this  case  one  has 

Theorem  2:  For  every  k,  there  exists  Ck  >  0  such  that,  if  f  1 , — is  a  multiplier  for 
the  fc-step  method  (p,  a )  then,  as 

m  >  c*n2-'. 

showing  the  blow-up  of  the  error  if  p  >  2  and  jj-»0. 

3)  A  graphical  method  was  devised  for  determining  whether  a  method 
possesses  a  multiplier  of  the  Popov  type.  This  type  of  multiplier  is  especially  conven¬ 
ient  for  proving  nonlinear  stability  results.  The  method  was  applied  to  exhibit  the 
stability  of  the  BDF  of  orders  2-5. 

C.  Nonlinearities  and  Multipliers 

Positivity  results  for  the  "energy-like"  quantities  of  the  form 

£  <  p.v„,  F„>  Z0,  (28) 

N 

where  Fn  ■  f(nh,  vn)  or  F„  «  f(nh,  x(nh)  +  v„)-f(nh,  x(nh ))  and  where  x(nh) 
denotes  the  exact  solution  of  the  stiff  system  and  <f>  the  scalar  product,  are  neces¬ 
sary  for  proving  nonlinear  stability  of  the  LMM  when  applied  to  x  •  f(t,  x).  We  give 
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examples  of  such  results  which  may  be  thought  of  as  "Generalized  Correlation 
Inequalities"  and  which  are  of  general  interest. 

Lemma  L 

1)  Let  d  be  convex,  non-negative  and  4(0)  m  0.  Let  (ij  £  0  for  1  and 
let  P  be  the  sequence  of  partial  sums  of  Ini;  then 

N 

2  <  #*•»'„.  grad  ♦„  >  a  (JW)*. 

2)  The  nonlinearity  is  said  to  be  gradient-like,  or  more  precisely, 
e-angle-bounded,  if 

</(*)  -  /O'),  y  -  z>  £  o<f(x)  -  /(r),  x  -  *>.  (29) 

An  example  are  the  3-monotone  functions  satisfying 

3 

2  <*/-*,- r  /(*,)>  £  0,x0**3- 

i-i 

In  this  case,  we  have 

Lemma  2:  If  /  is  o-angle-bounded,  /(0)  -  0  and  Hj  £  0  for  j  t  1 ,  then 

2  <**vn'fn>  *  +  ®>2  <V*A>-  (30) 

1 

3)  For  proving  convergence  results  for  LMM,  we  employ  the  following  useful 
inequality: 

Lemma  3:  Let  F(u,  v)  —  /(«,  x„  +  »)—/(«,  xH)  then 

2  F«>  *  M0-(i  Im,I  <v"’  f«>’ 
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where  Kv  K2  measure  the  asymmetry  and  time- variation  of  /,  respectively;  i.e.f 


<u,J(t)v>  4  -J-l <u,jy>  +  <>,»],  (31) 


<u,  J(l,  x)u>  4  K2  <«,  J(s,  y)u>,  (32) 


where  J  »  Jacobian  of  /. 

4)  Positivity  results  are  also  obtained  for  general,  not  necessarily  monotone, 
functions  if  they  are  nearly  linear,  as  in  the  following  lemma. 

Lemma  4:  Let  F(u,  v)  -  Hv  +  Hg(u,  v)  where  H  m  HT  2  0  and  assume  that 

<*(v),  Hg(v)>  4  L  <»,  Hv>. 


Then 


2  <*“V  F»>  2  <V  Hv„>, 

where  mQ  —  min  Re  p. 

D.  Boundedness  and  Stability  Results 

Various  boundedness  and  stability  results  for  the  numerical  solution  of 
difference  equations,  arising  when  A(a)-stable  methods  are  applied  to  nonlinear 
systems,  were  obtained  when  the  nonlinearity  and  the  multiplier  "match"  (in  the  sense 
that  the  inequality  (28)  holds,  as  is  the  case  for  example,  in  Lemmas  1-4  above). 

1.  Error  bounds. 
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We  consider  the  difference  equation 


p(£)x„  +  M£)F„  -  hpn+k  (33) 

where  pn  denotes  the  local  error  (or  a  forcing  term),  x„  the  global  error  (or  the 
computed  solution),  and  F„  the  global  error  in  /  (or  the  nonlinearity  itself,  respective¬ 
ly).  Then  if  F  satisfies  (28),  we  have  the  basic  boundedness  result  that,  for  some  C, 

|xj<  C[  max  |x  |  +/.[£,!  +  /»£  \p.\}.  (34) 

Equation  (34)  shows  that  global  errors  are,  under  such  conditions,  controlled  by  local 
errors. 

2.  Boundedness. 

We  have  shown  that  a)  if  the  perturbations  in  (33)  are  summable,  the?  (28) 
holds;  b)  if  /(0)  «  0  and  /  is  maximally  monotone  (though,  possibly,  discontinuous) 
then,  as  every  solution  of  the  difference  equation  tends  to  a  unique  value  of 

/“'(O).  Hence,  if  /"‘(O)  *  0  as  is  usually  the  case,  all  solutions  tend  to  zero  as  /-*■ ». 

3.  Convergence. 

It  was  shown  that  if  a  method  of  order  p  has  a  multiplier  compensating  for  the 
asymmetry  of  the  problem,  in  the  sense  that  |  pj  \  <  pQ  [see  equation  (31)],  then 
for  small  h  and  nh  <,  T  one  has 

|  global  numerical  error  |  £  CThp. 

4.  Input-output  stability. 

We  have  also  proved  f2-  and  fM-input  -  output  stability  of  (33)  under  slightly 


technical  assumptions. 


ffl.  EXTENDED  JOSEPHSON  JUNCTIONS  AND  INTERFEROMETERS 


A.  Background 

The  Josephson  phase  equation,  describing  an  extended  junction,  is 


,  d2*  d2*  d+ 

*  dx2~  9t2~  * 


k  sin 


(35) 


where  <f>  is  the  phase  and  0,  k  are,  resp.,  related  to  the  normal  and  Josephson  tunneling 


9&  5* 

currents.  One  seeks  solutions  4  such  that  —is  periodic  of  a  certain  period  where 


u  is  equal  to  the  d.c.  voltage.  If  these  solutions  are  non-static  they  are  called  running 
solutions.  The  boundary  conditions  on  <t>  depend  on  the  manner  in  which  the  junction 
is  driven  and  they  read 


♦  aJ1*-.-  ~H‘  <36> 

in  the  voltage-driven  case,  or  m 

I  *-0  -  -<*.  +  U  I  -  -  »e  (37) 

in  the  current-driven  case.  Here  He,  /,,  u  are  the  external  magnetic  field,  total 
current  and  d.c.  voltage,  respectively. 

A  "point”  junction  is  one  that  is  "small  enough"  for  the  term  92+/9x2  in  (35) 
to  be  negligible;  in  this  case,  the  phase  equation  becomes  an  ordinary  differential 
equation.  Interferometers  consist  of  a  given  number  n  of  point  junctions  connected  in 
various  manners  and  "loaded"  by  currents.  We  have  been  interested  in  studying  the 
properties,  especially  the  I-V  characteristics,  of  the  extended  junction  as  well  as  the 


symmetric,  two- junction  current  driven  interferometer,  whose  equations  are 


+  sin  +  A  m 

+  a#j  +  sin  $2  +  *  (^2“^i^ " ^/2  +  Ift 


(38) 


here  A  is  proportional  to  the  inductance,  Ic  is  a  control-current  and  <p  •  ut  +  periodic 
of  period  2 »/«,  where  u  is  proportional  to  the  voltage. 

B.  Extended  Junctions 

1.  Perturbation  results. 

The  running  solutions  and  I-V  characteristic  of  an  extended  Josephson  junc¬ 
tion  in  a  state  near  the  ohmic  regime  were  calculated  by  a  perturbation  method.  Both 
the  voltage-driven  and  current-driven  cases  were  considered  and  the  convergence  of 
the  perturbation  procedure  was  proved.  An  integral  representation  of  the  first  correc¬ 
tion,  in  the  I-V  curve,  to  the  ohmic  regime-as  well  as  its  dependence  on  the  external 
magnetic  field-was  given  and  evaluated  numerically  for  various  values  of  the  junction 
parameters.  The  details  of  these  results  were  given  in  [6]. 

2.  Existence,  uniqueness,  and  stability  results. 

If  one  writes  ♦  »  $0  +  where 

♦0  -  ut-kx  +  px2 


satisfies  the  formal  limit  equation  A$0  m  0  as  k-»0,  and  where  k  m  ou  +  Ht, 
p  m  ou/ 2,  and  +  is  periodic,  then  in  the  voltage  driven  case  +  satisfies 


A*  -«  sin  (♦„  +  *), 

(39) 

*  I  jr-0  “  ~jj~  l  "  °- 
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The  following  results  were  proved  in  this  case: 


a)  For  any  a«*0  there  exist  (2*/w-periodic  solutions  of  the  problem  (39). 
Hence,  for  or *0,  there  exist  (possibly  multiply  branched  I-V  characteristics  for  the 
Josephson  junction. 

b)  For  moderate  amounts  of  dissipation  relative  to  the  strength  of  the  nonlin¬ 
earity  (i.e.,  moderately  large  values  of  a  relative  to  a),  the  periodic  solution  of  (39)  is 
unique  and  globally  asymptotically  stable.  Hence  solutions  with  arbitrary  initial  data 

tend,  exponentially  in  t,  to  4  »  $0  +  t  with  ^  periodic.  For  example,  if  k  m  — 

8 

(corresponding  to  L/\j~. 9)  then  o>9/16  is  sufficient  for  uniqueness;  similarly,  with 
2  1 

k  m  1/4*  (corresponding  to  L/\j  m  uniqueness  is  assured  for  o>  .06. 

The  existence  theorem  a)  above  was  also  extended  to  the  current-driven 
junction  where  the  period  a  is  also  unknown. 

3.  Quasi-periodic  solutions. 

We  studied  the  behavior  of  the  extended  junction,  described  by  equation  (35), 
when  it  is  driven  by  an  oscillating  voltage  source.  More  specifically,  we  replaced  the 
first  of  the  boundary  conditions  (36)  by  the  condition 

♦  Ijc-o-PW  <40> 

where  p{t)  is  a  quasi-periodic  function  of  time  with  arbitrary  fundamental  frequencies 
«0,  <0,,...,  Such  problems  occur  if  the  junction  is  used  as  a  detection  device. 
Using  some  delicate  perturbation  and  estimation  arguments,  we  proved  that  for  o> 0 
the  phas^  equation  has  a  quasi-periodic  solution  which  is  unique  when  a  is  small. 

C.  Interferometers. 

1.  Existence  of  solutions 
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We  gave  a  proof  of  the  existence  of  running  solutions  to  (38)  for  arbitrary  k, 
under  the  assumption  that  the  dimension  of  the  connection  matrix  is  equal  to  one. 
The  proof  uses  a  homotopy  argument  to  show  that,  for  u#0,  oytO,  there  always  exists 
an  odd  number  of  solutions. 

2.  Numerical  perturbation  methods  for  interferometers 

We  developed  efficient  numerical  methods  for  finding  solutions  to  (38)  and 
the  associated  I-V  curves.  First,  for  a  given  /,  perturbation  methods  are  used  to  find  a 
phase  <f>  and  voltage  u  for  small  ac.  The  problem  is  then  formulated  as  a  boundary 
value  problem  for  the  pair  (4>,  u)  with  an  extra  condition  on  $  -  derived  from  the 
autonomous  character  of  (38)  which  "determines"  u.  Under  certain  assumptions,  we 
proved  that  We  pair  (<j>,  u)  is  isolated  and  that  a  Newton-like  method  converges.  The 
basic  numerical  procedure  is  a  continuation  method  with  gives  the  whole  I-V  curve. 
The  procedure  was  modified  in  order  to  improve  flexibility,  accuracy  and  efficiency. 
Specifically: 

a)  We  employed  a  two-parameter  continuation  method  to  compute  the 
solutions,  the  two  parameters  being  the  voltage  and  the  strength  of  the  nonlinearity. 
This  was  very  useful  for  finding  different  branches  of  solutions  by  following  different 
continuous  paths. 

b)  Since  the  calculation  was  delicate,  we  checked  its  accuracy  by  three 
different  methods: 

1 )  by  grid  refinement; 

2)  by  deriving  a  consistency  relation  /  ■*  I(4>)  which  the  exact  solu¬ 
tion  must  satisfy  and  monitoring  the  discretization  error  in  this  relation;  and 

3)  by  carrying  out  linear  stability  calculations. 

(Since  the  translation  invariance  of  equation  (38)  implies  that  at  least  one  Floquet 
multiplier  must  be  unity,  this  fact  was  also  used  to  check  the  accuracy.) 
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c)  To  obtain  a  complete  picture  of  the  interferometer,  the  basic  calculation 
has  to  be  repeated  104~105  times  and  thus  must  be  done  quite  efficiently.  The 
following  are  some  of  the  features  which  enhanced  the  efficiency  of  our  computation: 

1)  The  unknowns  in  the  linearized  equations  are  numbered  in  an 
appropriate  way  for  sparse  matrix  techniques  to  become  applicable.  Since 
these  equations  depend  strongly  on  «,  and  weakly  on  X,  the  computation  of 
"X— w  grid"  was  arranged  so  that  the  outer-most  loop  was  the  w-loop,  the  next 
was  the  X-loop  and  the  inner-most  was  the  Newton  loop. 

2)  The  computation  took  advantage  of  the  consistency,  in  X,  of  the 
off-diagonal  part  of  the  linearized  equations. 

3)  The  starting  guess  for  the  Newton-loop  was  obtained  by  using  a 
"predictor";  this  reduced  the  number  of  the  Newton-steps  needed  for  conver¬ 
gence. 

3.  Numerical  results. 

The  above  numerical  study  yielded  the  following: 

a)  The  correct  dependence  of  the  resonant  current  on  the  two  physically 

important  parameters  X  and  Ta—^* 

V2 

physical  analysis  gave  the  dependence  of  this  resonant  current  on  T  alone.  Our 
calculations  showed  that  these  previous  results  are  qualitatively  valid  only  for  X  around 
0.8- 1.0.  However,  they  underestimate  the  maximum  resonant  current  by  about  SO 
percent  for  these  values  of  X  and  are  rather  invalid  elsewhere. 

b)  There  are  subharmonic  resonances  in  the  I-V  curve  which  occur  around 
one-third  of  the  resonant  frequency  (voltage)  for  small  X,  but  shift  towards  one-half  of 
the  resonant  frequency  for  X~l,  i.e.,  in  the  strongly  nonlinear  case. 

c)  At  resonance,  the  two  phase  functions  are  opposite  in  sign  and  for  each 
function  all  the  Fourier  components  are  negligible  except  for  the  fundamental  frequen- 


Previous  approximate  calculations  and 
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cy.  However,  for  small  voltage,  ~.l  of  resonance,  the  phase  functions  have  an  almost 
saw-tooth  shape  with  many  frequencies  present. 

d)  Two  types  of  bifurcations  for  the  solutions  to  (38).  Near,  and  a  little  to 
the  left  of  the  main  resonance,  bifurcation  to  another  (2tr)u)-periodic  solution  occurs. 
Bifurcations  to  (4ir)«)-periodic  solutions  occur  to  the  right  of  the  main  resonance  (for 
values  of  A~2). 

e)  The  dependence  of  the  amplitude  of  the  subharmonic  resonances,  as  well 
as  the  main  one,  on  the  strength  of  the  nonlinearity  (i.e.,  on  the  parameter  X  in  (38)), 
was  computed.  The  subharmonics  are  about  IS  percent  of  the  main  resonance;  they 
are  stable  for  small  nonlinearities  (X~.S)  but  unstable  at  A~l. 

f)  The  1-V  characteristic  curves  were  found  to  be  generally  multi-valued  and 
to  possess  many  loops  (even  though  most  of  the  loops  are  either  unstable  or  only 
meta-stable).  Hysteretic  behavior  is  clearly  evident. 

g)  Negative  resistance  portions  of  the  I-V  curve  occur  near  the  main  reso¬ 
nance  \=«l).  These  portions  are  stable,  perhaps  contrarily  to  physical  intuition. 

We  also  calculated  the  (linear)  stability-instability  behavior  of  all  of  the 
computed  branches  of  solutions  to  (38). 
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IV.  SPECIAL  METHODS  FOR  QUADRATIC  SYSTEMS  OF  DIFFERENTIAL 
EQUATIONS 

The  fractional  linear  scheme 

u(:  +  h)  m  (-/?(*)[«</)  1  +  S(*))''(T(/i)u(t)  +  v(h ))  (41) 

was  previously  applied  to  the  numerical  solution  of 

i  -  q{tjc) 

where  q  is  quadratic.  It  was  shown  [11]  that  convergent  second-order  accurate 
schemes  of  the  above  form  exists  for  all  quadratic  systems  and  an  explicit  construction 
of  such  a  scheme  was  given. 

These  results  were  extended  to  equations  of  the  form 

x  =  (  -  KW  +  I)~\x»x  +  Bx  +  c)  (42) 

where  R[x]x  m  0.  In  particular,  it  was  shown  that  i)  a  convergent  fractional  linear 
scheme  converges  to  an  equation  of  type  (42);  ii)  every  equation  of  type  (42)  can  be 
solved  numerically  by  a  second-order  accurate  convergent  fractional  linear  scheme. 

A  complete  set  of  invariants  was  derived  for  all  two-by-two  systems  of 
first-order  quadratic  differential  equations,  with  respect  to  the  six-parameter  group  of 
all  affine  transformations.  These  invariants  provide  a  criterion  for  deciding  to  which 
affine  equivalence  class  any  given  system  belongs.  Furthermore,  all  menbers  of  any 
one  of  these  equivalence  classes  were  shown  to  be  equivalent  to  a  single  second-order 
equation  in  one  dependent  variable.  Special  cases  of  the  latter  define  the  set  of  all 
elliptic  functions  and  other  abelian  integrals.  A  novel  convergent  second-order 
discretization  of  the  above  mentioned  second-order  equation  was  also  proposed. 
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